#include <cmath>
#include <cstdio>
#include <iostream>
#include <vector>
#include "matrix.h"
#include "MatrixRegression.h"  

#ifndef _NO_NAMESPACE
using namespace std;
using namespace math;
#define STD std
#else
#define STD
#endif

#ifndef _NO_TEMPLATE
typedef matrix<double> Matrix;
#else
typedef matrix Matrix;
#endif

#ifndef _NO_EXCEPTION
#  define TRYBEGIN()	try {
#  define CATCHERROR()	} catch (const STD::exception& e) { \
						cerr << "Error: " << e.what() << endl; }
#else
#  define TRYBEGIN()
#  define CATCHERROR()
#endif


double Regression
(double underly, vector<double> Strke, vector<double> Time, vector<double> vols, vector<double>& Bs)
{
	TRYBEGIN()
	{
	    int vectorsize = Strke.size()-1;

	    //Convert our vectors into Matrix form.
		Matrix X(vectorsize, 3);
		Matrix Y(vectorsize, 1);
		Matrix betas(1, 5);
		vector<double> vec_of_betas(3);
			
		int choice;
		double B0, B1, B2, B3, B4;
		double S,K,T,T2;
		double vol_estimate, vol_estimate2, var_estimate, var_estimate2, forward_vol;
		
        for (int i = 0; i < vectorsize; i++) {
            X(i,0) = 1;
			X(i,1) = log(Strke[i]/underly)/sqrt(Time[i]/365);
            X(i,2) = pow(log(Strke[i]/underly)/sqrt(Time[i]/365),2);
            Y(i,0) = vols[i];
        }
                 
        // computes coefficients using Least Squares Approximation
        betas = !(~X*X)*~X*Y;
        B0 = betas(0,0);
        B1 = betas(1,0);
        B2 = betas(2,0);

	  Bs.push_back(B0);
	  Bs.push_back(B1);
	  Bs.push_back(B2);
        
        cout << "Beta0 = " << B0 << endl;
        cout << "Beta1 = " << B1 << endl;
        cout << "Beta2 = " << B2 << endl;
        
        cout << "The equation to approximate sigma is: " << endl;
        cout << "sigma = " << B0 << " + " << B1
			 << "*log(K/S)/sqrt(T) + "<<B2<<"*[log(K/S)/sqrt(T)]^2" << endl;
        
		// If you want test the functions, just delete the //
        // printDebug(B0, B1, B2); 
	}

	CATCHERROR();
}

double getVolatility(double k, double s, double t, double beta0, double beta1, double beta2)
{
	double vol_estimate = beta0 + beta1*(log(k/s)/sqrt(t/365))
							    + beta2*(log(k/s)/sqrt(t/365))*(log(k/s)/sqrt(t/365));
	return vol_estimate;
}

double getVariance(double k, double s, double t, double beta0, double beta1, double beta2)
{
	double vol_estimate = beta0 + beta1*(log(k/s)/sqrt(t/365))
							    + beta2*(log(k/s)/sqrt(t/365))*(log(k/s)/sqrt(t/365));
	double var_estimate = vol_estimate * vol_estimate;

	return var_estimate;
}

double getForwardVolatility
(double k, double s, double t, double T, double beta0, double beta1, double beta2)
{
	double vol_estimate = beta0 + beta1*(log(k/s)/sqrt(t/365))
							    + beta2*(log(k/s)/sqrt(t/365))*(log(k/s)/sqrt(t/365));
	double vol_estimate2 = beta0 + beta1*(log(k/s)/sqrt(T/365))
							    + beta2*(log(k/s)/sqrt(T/365))*(log(k/s)/sqrt(T/365));
	double var_estimate = vol_estimate * vol_estimate;
	double var_estimate2 = vol_estimate2 * vol_estimate2;
	double forward_vol = sqrt((var_estimate2*T - var_estimate*t)/(T-t));

	return forward_vol;
}

void printDebug(double beta0, double beta1, double beta2)
{
	cout << "This is a test of calling getVolatility: "
		 << getVolatility(100, 110, 200, beta0, beta1, beta2) << endl; 

	cout << "This is a test of calling getVariance: "
		 << getVariance(100, 110, 200, beta0, beta1, beta2) << endl; 

	cout << "This is a test of calling getForwardVolatility: "
		 << getForwardVolatility(100, 110, 200, 300, beta0, beta1, beta2) << endl; 
}
